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(N ■ Abstract 

We calculate the chiral condensate of QCD with 2, 2+1 and 3 flavors of sea quarks. Lattice 
QCD simulations are performed employing dynamical overlap fermions with up and down quark 
^ ' masses covering a range between 3 and 100 MeV. On L 1.8-1.9 fm lattices at a lattice spacing 

~ 0.11 fm, we calculate the eigenvalue spectrum of the overlap-Dirac operator. By matching the 
lattice data with the analytical prediction from chiral perturbation theory at the next-to-leading 
order, the chiral condensate in the massless limit of up and down quarks is determined. 
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I. INTRODUCTION 



Chiral condensate E in Quantum Chromodynamics (QCD) is not a directly accessible 
quantity in experiment, yet plays a crucial role in the low-energy dynamics of QCD as 
an order parameter of chiral symmetry breaking. When E in the limit of massless quarks 
is nonzero, chiral symmetry is spontaneously broken and hadrons acquire a mass of order 
^QCD; the QCD scale. Only the pion remains massless as the Nambu-Goldstone boson; 
its dynamics is well described by an effective theory known as chiral perturbation theory 
(ChPT) [l, 2|. E is one of two most fundamental parameters in ChPT (the other is the pion 
decay constant F) and appears in the predictions of various physical quantities. 

Calculation of E, the expectation value of the scalar density operator E = —{qq), from 
the first principles of QCD requires nonperturbative techniques. In this paper we report on 
a numerical calculation of E in lattice QCD including the effects of up, down and strange 
sea quarks. We investigate the low-lying eigenvalue spectrum p(A) of the Dirac operator, 
which is related to E through the Banks-Casher relation p(0) — )■ E/vr jsl and its extension 
to the case of nonzero A. Since only the low-lying eigenvalues are relevant, one can avoid 
contamination from ultraviolet divergence of the scalar density operator qq, which is of order 
m/a^ at a finite quark mass m and a lattice spacing a. 

The Banks-Casher relation is satisfied in the limit of massless quarks after taking the 
large volume limit (the thermodynamical limit), which is the meaning of the arrow in the 
relation p(0) — )■ E/vr. When the massless limit is taken at a finite volume, the vacuum 
expectation value of qq shows a critical fluctuation which leads to a divergent correlation 
length and vanishing chiral condensate. Taking the thermodynamical limit, on the other 
hand, is numerically expensive in practical lattice calculations. In this study, we use the 
low energy effective theory as a guidance of volume and quark mass scalings. Namely, once 
the scaling behavior predicted by the effective theory is confirmed by the lattice data, the 
infinite volume and chiral limit according to the scaling can be safely taken. The scaling we 
consider is that for varying the eigenvalue A, volume V and the quark mass m. 

We use the ChPT formula for the low-lying eigenvalue spectrum of the Dirac operator 
{4], which is valid in both the p and e regimes. The e regime is a parameter region where the 
quark mass is so small that the Compton wavelength of the pion is longer than the extent of 
the finite-volume space-time. In this regime, the constant mode of the pion field has to be 
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9|, llO| . The new method 



integrated out over the group manifold when the path integral is evaluated, which is a non- 
perturbative calculation in the sense that one does not use the expansion in small pion field. 
At the leading order of the so-called e expansion, the system is equivalently described by the 
chiral Random Matrix Theory (ChRMT), with which a number of theoretical predictions 
for the low- lying Dirac spectrum have been derived 5|-|8|. In the other regime (p regime), 
where the pion Compton wavelength fits in the volume, the conventional ChPT applies and 
the calculation of the Dirac operator spectrum is available as well 

given in j4| consistently combines the both results within a systematic expansion, and thus 
is valid in both regimes as well as in between. By lattice calculations, we produce the data 
at various sets of quark masses to firmly test this analytic expectation. 

Application of the ChPT formula for the low-lying Dirac spectrum requires a good control 
of the chiral symmetry in the calculation. In this work, we use the overlap-Dirac operator 
Q Q wh,ch satisfies the G,„sparg-Wilso„ relation Q and thus .eal^es a .noMfie, chiral 
symmetry on the lattice [M]. Since the chiral effective theory is constructed only assuming 
the presence of chiral symmetry, the same construction as in the continuum theory can be 
applied to lattice QCD with overlap fermions. Although the overlap-Dirac eigenvalues lie 
on a circle on the complex plane, the physical imaginary part is uniquely identified up to 
0{a?) discretization effects. 
In our previous works 
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16| we performed large-scale lattice simulations of two-fiavor 



QCD using the overlap fermion formulation 17|] and calculated low-lying Dirac eigenvalues. 
By matching the lowest eigenvalue spectrum with the ChRMT expectations, we extracted 
S. Since the ChRMT corresponds to the leading order of the e expansion in ChPT, this 
result is subject to NLO or 0{1/ F'^V^^'^) corrections, which could be sizable on the lattice 
of size V = L^T ~ (1.9 fm)^ used in that study. Another limitation was that the quark 
mass must be very small to apply the e regime formula, and the runs in the p regime could 
not be used in the analysis. 

Application of the new formula jj] is attempted for the first time in our recent work jisl-fioll 
in which next-to-leading order (NLO) corrections are included in the analysis. Using the 
2+1-fiavor QCD data generated with dynamical overlap fermions on 16'^ x 48 lattices, the 
value of S in the chiral limit of two light quarks is obtained. The present paper provides an 
extensive description of this work. 

In this paper, we analyze the low-lying Dirac spectrum mainly on the 2-|-l-fiavor QCD 
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simulations, where the strange quark mass is fixed near its physical value. The mass 
of degenerate up and down quarks, rriud, covers a range between m^/S and nis in the p 
regime lattice ensembles, which are generated for calculations of various physical quantities 
including the pion mass and decay constant 20|. We also generate an e regime ensemble, 
where the up and down quarks are kept nearly massless while is fixed near the physical 
value. The NLO ChPT formula allows us to combine these data to obtain S in the chiral 
limit of up and down quarks. We also extract the pion decay constant F and one of the 
NLO low energy constants (LECs) Lq from the correction terms. 

As demonstrated in the following sections, the ChPT formula provides information on 
the shape of the low-lying Dirac spectrum, with which we can test the agreement between 
the formula and the lattice data in detail. The volume dependence gives a critical test, since 
it is essentially controlled by S and F. At some parameter points, we compare the data 
obtained on a larger (24^ x 48) lattice to those on the smaller (16^ x 48) volume to check if the 
ChPT prediction consistently describes the difference. The sea quark mass dependence of 
the chiral condensate is partly controlled by Lq but we also expect a nonanalytic dependence 
due to the pion-loop effects that should be observed in the lattice data. These nontrivial 
consistency checks are performed to gain confidence in the final result for S. 

In addition to the main analysis in 2+1-fiavor QCD, we also investigate two-fiavor QCD 
and three-fiavor QCD where rriud and are degenerate. For the case of two-fiavor QCD, 
the lattice data are the same as in our previous studies 16], but we use the ChPT 
formula valid at NLO in this new analysis. The degenerate three-fiavor QCD configurations 
are newly generated for this study at two light quark masses. We thus obtain the chiral 
condensate S for these variants of QCD. 

This paper is organized as follows. In Section [Tll we briefiy explain how the chiral con- 
densate is determined from the Dirac eigenvalue density and how the finite volume effects 
are removed using ChPT. Some technical aspects of lattice QCD simulations are described 
in Section UTTl and the numerical results are discussed in the following sections. First, we de- 
scribe the strategy to extract the LECs from lattice data of the spectral density in Section HVl 
Second, numerical scaling tests of the NLO ChPT formula are presented in Section |Vl We 
then proceed to the determination of S in the chiral limit in Section IVI[ Summary of this 
analysis and conclusions are given in Section IVIIi 
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II. BANKS-CASHER RELATION IN A FINITE VOLUME 



Eigenvalues of the Dirac operator in four- dimensional continuum Euclidean space are 
pure imaginary (which we denote iXi, 1X2, ■ ■ ■ with real Aj's). The spectral density at an 
eigenvalue iA is defined by 

-| 00 

^ tfEW-^^))' (1) 

k 

where (■ ■ ■) denotes an average over the gauge configuration space. Since a nonzero eigen- 
value appears as a pair with its complex conjugate, p(A) = p(— A), we only consider the 
A > region in the following. 

The chiral condensate S in the limit of massless quarks and infinite volume is an order 
parameter of the chiral symmetry breaking. Through the Banks-Casher relation [3], S is 
related to p(A) as 

lim lim p(0) = -, (2) 

with which one can identify the spontaneous breaking of chiral symmetry by measuring p(0) 
instead of S. 

Even when the volume V, sea quark masses {m^ea} = {^u,'^d,''TT's, ' ' '} ^ind A are all 
finite, a similar nonperturbative relation 



p(A) = -Re- 



71 



(3) 

m,y=iX 



holds. Here, (■ ■ ■)|m.„=jA means that the quantity is evaluated with the valence quark mass 
analytically continued to a pure imaginary value iX. In this relation, the ultraviolet diver- 
gence in the definition of the qq operator cancels by taking its real part at an imaginary 
value = iX (where the divergent part is pure imaginary), which is natural because the left 
hand side of the equation only refers to low-lying modes and is insensitive to the ultraviolet 
region of the dynamics. 

We note that the relation ([3]) is valid even when the ensemble average (■ ■ ■) is restricted 



to a given topological sector of Q in the gauge field configurations [2l|], that we denote 
(■ ■ ■)q. Namely, if we define the spectral density at a given topological sector as Pq(A) = 
(l/V) J2k {^{^ '~ ^k))Q, it is obtained by computing —Re{qq)Q/7i. 

Although the chiral condensate E is different from — Re(gg)|m„=iA at finite volumes, the 
difference can be described by the low-energy effective theory, provided that the energy 



scales of the theory are well below the QCD scale: 

A, m„ XIV^I^ « Aqcd. (4) 

It is, therefore, possible to directly compare the lattice QCD calculation of p(A) at finite 
A, mj, Q and V with the prediction of the effective theory. By taking the limit of mj — )■ 
after V ^ oo according to the effective theory and summing over Q, one can reproduce the 
physical p(A), which in the limit of A — t- gives E. 



I. th,s direction, studies in both latfce QCD the low-energy effective theory 

have been done. Smilga and Stern [9|] and Osborn et al. [10| calculated the Dirac eigenvalue 
spectrum in the conventional p expansion of (partially quenched) ChPT to NLO. In the 
vicinity of A = 0, which corresponds to the limit of zero valence pion mass, a special 
treatment of the zero-momentum modes is needed because the correlation length exceeds 
the size of the volume. This special treatment is known as the e expansion of ChPT, in 
which the zero-momentum modes are nonperturbatively integrated over the SU{Nf) (or 
U{Nf) when the topological charge Q is fixed) manifold. At the leading order (LO), the 
finite size effect around A ~ was calculated in . Their results are expressed using the 
Bessel functions, which has a ~ gap from zero, reflecting the fact that no spontaneous 

symmetry breaking occurs at finite volumes. ^ 

Recently, an interpolation between the p and e regimes was considered in [4!| . The recipe 
for the calculation is to keep the same counting rule as in the p expansion but to integrate 
the zero-modes nonperturbatively like in the e expansion. Partial quenching is performed 
with the so-called replica trick so that results at arbitrary nondegenerate set of quark masses 
can be compared to lattice QCD. Using this hybrid method, the results mentioned above (in 
the p expansion j^, [l^ and in the e expansion 0-^) are smoothly connected. Comparison 
with the lattice data is, therefore, no longer limited in either the e or p regimes, and more 
precise determination of S is possible. 

Here we briefly reproduce the result of 4] where we consider a general theory with Nf 
flavors of sea quarks. The spectral density in a fixed topological sector of Q is given by 

Pq(A) = ScSPQ(AScffl^, {nisea^esV}) + ff{\, {nisea}), (5) 

where two terms pq^XTiesV, {rrisea^esV}) and //'(A, {nisea}) are given in the following. Seff 
includes the leading finite quark mass correction to S that modifies the overall normalization 
of the spectrum, and is therefore called the effective chiral condensate. 



The spectrum of the near- zero modes (A ~ l/HV) is mainly affected by the zero- 
momentum pion fields. The first term in ([5]) has the same functional form as the one 



at the leading order of the e expansion 



combinations XTj^sV and {nisea^csy} = {'^iSgfrV, 

ICI 



,(C,{/^sea}) = C2- 



which is expressed in terms of dimensionless 
detB 



(6) 



with Nf X Nf matrix A and {Nf + 2) x [Nf + 2) matrix B defined by Aij = fi{ ^lQ+j^i{fii 



and B,, = C^-Vq+,_2(C), = C^--Vq+,_i(C), 4 



-f^t-2y (i 7^ 1,2), 



respectively (J^'s and J/'s denote the (modified) Bessel functions.). The phase factor C2 is 
1 for Nf = 2 and 3. 

The second term in ([5]) is a logarithmic NLO correction (chiral-logarithms) which is also 
partly seen in the conventional p expansion [10]. With the meson mass M?- = (mj+mj)E/F^, 
which is made of either sea quark (/) or valence quark (v), it is given by^ 



ff{X, {rUsea}) 

where 

1 



Re 



Nf 



E(A(M|J - A(Mj^/2)) - (G(M- 
/ 



vv) 



am 



(7) 



A(M2) + - M2,)9m2A(M2 



(Nf = 2), 



A(M2 



< 1 

3 



M2 
16^ 



2(M, 



ud 



2 

ss) 



9(M2 - M2)2 

, (M2-M2,)(M2 



(M2 - M2) 



9(M2 - M2)2 ^ 

M2 1 

-^9m2A(M2) 



A(M2) 



(8) 



1, 3), 



In 



M2 



+ ^i(M2). 



(9) 



Here, the physical meson masses are given by the leading order relations M^^ 



2m„S/F2 = 

2mrfS/F2, Ml = 2m,S/F2 and = {M^^ + 2MfJ/3. The scale fXsub (= 770 MeV in this 
work) is a subtraction scale. The function given by gi{M'^) = gi{M'^) — 1/M21^ denotes 
a finite volume correction from nonzero momentum pion modes. In the p expansion, it is 
expressed by the modified Bessel function Ki 29| while in the e expansion a polynomial 



^ In this paper, we use simplified notations: A(M^) and G(Af ^) correspond to A(0, KP) and 0(0, hP, M^) 
in respectively. 
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expression is used [30|. In this study, we need the both expressions: 



E 



VM2 
47r2|a| 



Ki{VM^\a\) - 



1 



{\M\L > 2) 



(10) 



ln(MV^/2) _ J2 



n=l 



(n 



1)! 



j\^2(n-l)p.(n-2)/2 < 2) 



where denotes a four- vector 



shape coefficients defined in |3(j]. Their formula and numerical values for the first several /3„'s 
are summarized in Appendix |X1 In our numerical study, we truncate the sum at n^"-^ = 7 



[niL,n2L,n3L,n4T) with integer nj's and /3j's are the 



and n^'^^ 



300, which indeed shows a good convergence around the threshold \M\L = 2. 
We note that both A(M^) and G'(M^) are finite even in the limit of M — )■ 0. 
The effective chiral condensate Sgfr in ([5]) is given by 



^cff 



'Nf 



Nf 



1 - I E A(Mj^/2) - G(0) - 16LI Y: Mj 



(11) 

/ / ' ' 

where Lg (renormalized at figub) is one of the low-energy constants at NLO 2|. From the 
sea quark mass dependence of Seg, one can determine S as well as F and Lg. 

In the expression (jS]), dependence on the topological charge Q is encoded only in the first 
term T^esPgi^^csV, {rrisea^csV}) and the second term ff{\, {rrisea}) does not depend on Q 
since it is a contribution from nonzero momentum modes. On the other hand, the chiral 
logarithm manifests itself in the both terms through A(M^). Since could also contain 
A through m„ = iX, the spectral density shows a nonanalytic functional form. 

The above ChPT results are subject to higher order corrections in the p expansion, 
for which the expansion parameter is either M'^/F'^ or l/(FV^^^Y. Although the zero- 
mode contribution is treated nonperturbatively, there are two-loop contribution of nonzero 
momentum modes that could also couple to the zero-mode and introduce different types 
of group integrals. At the two-loop level, these contributions may have the order M^/F^, 
M"^ / [F^V^^"^) or 1/{F^V), whose coefficients are unknown. We therefore need to carefully 
check the convergence of the expansion at NLO for our parameter sets. In the following 
analysis, we test the NLO formula with various sets of quark masses, as well as different 
lattice volumes {L = 16, 24) for the Nf = 2 -|- 1 runs and different topological sectors for 
the Nf = 2 runs, in order to confirm the convergence. 
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III. LATTICE QCD SIMULATIONS 

Numerical simulations of lattice QCD are performed with the Iwasaki gauge action [3V\ 
at /3 = 2.3 (except for the run of Nf = 2 QCD at m^^ = 0.002 for which we choose /3 = 2.35) 
including 2, 2+1 (m^ fixed), and 3 (m^^ = m^jflavors of dynamical quarks. For the quark 



action, we employ the overlap-Dirac operator 11 | 



D{m) = (mo + Y ) + ( "^0 - y j 75Sgn[ifH/(-mo)], (12) 

where m denotes the quark mass and ifvK = 75-Dh/(— ^o) is the Hermitian Wilson-Dirac 
operator with a large negative mass —mo. We take mo = 1.6 throughout our simulations 
(here and in the following the parameters are given in the lattice unit.). For the details of 
numerical implementation of the overlap-Dirac operator, we refer to our previous paper 

It is known that the numerical cost for the dynamical simulation of the overlap fermions 
becomes prohibitively large when Hw{—mQ) has (near) zero- modes. To avoid this problem, 
we introduce extra Wilson fermions and associated twisted mass bosonic spinors to generate 
a weight 



det[g^(-mo)] 



(13) 



in the functional integrals |2ll . 1321 . |33|. Both of these fermions and ghosts are unphysical 
as their masses are of order of the lattice cutoff, and do not affect low-energy physics. The 
numerator suppresses the appearance of near-zero modes, while the denominator cancels 
unwanted effects from high modes. The twisted-mass parameter mt controls the value of 
threshold below which the eigenmodes are suppressed. In our numerical studies, we set mt 
= 0.2. 



With the determinant (11311 the index of the overlap-Dirac operator, or the topological 



charge in the continuum limit |34| , never changes from its initial value during the molecular 
dynamics steps since its change always requires crossing zero eigenvalue of Hw{—mo). In 
this work the simulations are mainly performed in the trivial topological sector, Q = 0. In 
order to check the topological charge dependence, we also carry out independent simulations 
at Q = +1, —2 and —4 at several sets of parameters. 

Simulation parameters are summarized in Table [B The lattice size is = 16'^ x 32 for 
the Nf = 2 runs, while it is \^ = 16^ x 48 for the main Nf = 2 + 1 and Nf = 3 runs. In 
order to investigate the finite volume scaling, we also simulate on a 24^ x 48 lattice at the 
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Nf V P a"^(GcV) m„d Q Ntrj nrj Nauto 

2 16^ X 32 2.35 1.776(38) 0.002 oo 4680 0.5 34(12) 

2.30 1.667(17) 0.015 oo 10000 0.5 48(21) 

0.025 oo 10000 0.5 38(16) 

0.035 oo 10000 0.5 28(12) 

0.050 oo 10000 0.5 24(9) 

0.050 oo -2 5000 0.5 50(24) 

0.050 oo -4 5000 0.5 34(16) 

0.070 oo 10000 0.5 23(10) 

0.100 oo 10000 0.5 9(3) 

2+1 16^ X 48 2.30 1.759(10) 0.002 0.080 5000 0.5 17(9) 

0.015 0.080 2500 1.0 15(6) 

0.015 0.080 1 1800 1.0 5(1) 

0.025 0.080 2500 1.0 11(5) 

0.035 0.080 2500 1.0 24(11) 

0.050 0.080 2500 1.0 9(5) 



0.015 0.100 2500 1.0 5(3) 

0.025 0.100 2500 1.0 10(3) 

0.035 0.100 2500 1.0 34(21) 

0.050 0.100 2500 1.0 5(3) 



24^ X 48 2.30 1.759(10) 0.015 0.080 2500 1.0 2(1) 

0.025 0.080 2500 1.0 3(2) 

3 16^ X 48 2.30 1.759(10) 0.025 0.025 2500 1.0 4(1) 

0.035 0.035 2500 1.0 5(1) 

0.080 0.080 2500 1.0 20(12) 

0.100 0.100 2500 1.0 24(15) 

TABLE I: Summary of lattice parameters. For the 2-, 2+1- and 3-flavor runs, the values of 
/3, a~^, {i^udii^s) and Q are listed. Other parameters are those for the HMC simulations: Ntrj 
denotes the number of trajectory, Ttrj is the unit trajectory length, and Nauto denotes the integrated 
auto-correlation length (number of trajectories) of the lowest eigenvalue. 
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same lattice spacing with two choices of sea quark masses {mud^ms) = (0.015,0.080) and 
(0.025,0.080). For the determination of the lattice spacing a = 0.11-0.12 fm (a^^ = 1.7 - 
...8 GeV), we choose the f2-baryon mass as the input for the Nf = 2 + 1 and 3 ensembles 
351], while it is determined from the heavy quark potential with an input Vq = 0.49 fm for 
the Nf = 2 case. Our lattice size is then estimated as L ~ 1.8 fm for the Nf = 2 + 1 runs, 
and L ~ 1.9 fm for the Nf = 2 runs. 

In the Nf = 2 runs, seven different values of the up and down quark mass rriud are taken. 
For the Nf = 2 + 1 runs, we choose two different values of strange quark mass (= 0.080 
and 0.100) and six (for the former) or five (for the latter) values of rriud are chosen. For the 
degenerate Nf = 3 flavor runs, we take four different values of rriud = ttLs- Note that the 
lightest up and down quark mass rriud = 0.002 in the Nf = 2 or Nf = 2 + 1 runs roughly 
corresponds to 3 MeV in the physical unit, with which pions are in the e regime while kaons 
still remain in the p regime. 

We compute 50-80 pairs of low- lying eigenvalues (that we denote A°^'s) of the massless 
overlap-Dirac operator -D(O) at every 5 or 10 (depending on the parameters) trajectories. 
In the calculation of the eigenvalues, we employ the implicitly restarted Lanczos algorithm 
for the chirally projected operator P+ -D(O) P+ (where P+ = (1 + 75)/2) of which eigenvalue 
corresponds to ReA°''. From each eigenvalue of P+D{0) P+, the eigenvalue A''^, as well as 
its complex conjugate, are extracted through the relation |1 — A°"/moP = 1. In order to 
compare with ChPT, every complex eigenvalue A°^ is mapped onto the imaginary axis as 
A = ImA°''/(l — ReA°^/(2mo)). The difference between A and ImA"^ is a discretization effect, 
which is negligible (within 1%) for |A"''| < 0.03. In the analysis, we consider positive A only. 

For each run, 1,800-10,000 (depending on the parameters) trajectories are accumulated 
using the hybrid Monte Carlo algorithm. The integrated auto-correlation time Nauto of the 
lowest A in the unit of the trajectory length Tfrj is also listed in Table [H Because of its 
infra-red nature, the lowest Dirac eigenvalue is expected to be most difficult to decorrelate 
and thus has the longest auto-correlation time. The measurement is not stable and the 
statistical error is as large as 50% in some cases, but Nauto is typically 0(50) or less. In 
the following analysis, the statistical error for the spectral density and other quantities is 
estimated by the jackknife method after binning the data in every 100 trajectories. 

Details of configuration generation and other quantities will be reported in a separate 
paper. 



11 



IV. EXTRACTION OF LEGS AT EACH SET OF SEA QUARK MASSES 



Although a global fit of the lattice data for spectral density to the NLO ChPT formula is 
possible in principle, we prefer to simplify the analysis, for better understanding of numerical 
sensitivity of the lattice data and the errors in the final result. We first consider the mode 
number below A, or the integrated eigenvalue density, 



at each set of sea quark masses. The analytic ChPT result (jS]) is also integrated numerically 
from to A. In the second term of ([5]) we can replace S/F^ by Scff/-F^ as their difference is 
a higher order effect. Then, there are two unknown parameters in the formula: Seg and F. 

The data points of Nq{X) at two reference values of A are hence sufficient to determine 
the parameters. As the reference points we take A = 0.004 (~ 7 MeV) and 0.017 (~ 30 MeV) 
except for the case with rriud = 0.002, for which we choose A = 0.0125 and 0.017 {Nf = 2) 
or 0.010 and 0.017 {Nf = 2 + 1). For the Q 7^ runs we take A = 0.01 and 0.02. Effectively, 
the lower A point determines Scg, while the other point is more sensitive to the NLO effects 
that contain l/F"^. We check that the resulting values of Scg and F are stable against the 
change of the reference points by varying them by a factor of 2 or 3 while keeping the higher 
point less than 0.030 to avoid possible higher order corrections. The reference points are 
different for the e regime runs and for the Q 7^ runs, because the small eigenvalues are 
highly suppressed (the lowest eigenvalue is larger than 0.004) for these cases. 

For the Nf = 2+1 lattice data, we test both the Nf = 2 + 1 and Nf = 2 ChPT formulas. 
For the latter case, the strange quark is assumed to be decoupled from the theory, which we 
call reduced Nf = 2 ChPT. 

Numerical results are listed in Table [Til Before moving to further analysis of the results 
let us describe our observations for the spectral function. 

Figures [TH8] show the spectral density and its integral f|T^ obtained in our lattice sim- 
ulations. A typical example is that for Nf = 2 + 1 QCD in the p regime (Figured]). The 
upper panel shows the histogram plot of the spectral function Pq(A) as a function of A. The 
lattice data for each bin have a jackknife estimate of the statistical error. The solid (red) 
curve represents the NLO ChPT formula, while the dotted (blue) curve corresponds to the 
leading order result (in the e expansion). Since the first reference point is 0.004, it probes 




(14) 
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the first peak, which corresponds to the lowest eigenvalue. Starting from the second peak, 
the effect of the NLO term appears, as clearly seen from the difference of the two curves. 
Therefore, by taking the second reference point at 0.017, the data have enough sensitivity 
to the NLO parameter l/F"^. This observation is of course specific to the particular volume 
of our lattice; on larger lattices, the peaks move toward the origin and the impact of the 
NLO term would become less significant on the second or third peaks (see Figure [8]). 

The agreement with the formula can be seen more clearly by looking at Nq{X), the mode 
number below A (lower panel). The lattice data depart from the leading-order curve (dotted), 
which corresponds to the first term J^^sPq in dS]), at around 0.005 (see the inset). Then, the 
data follow the nontrivial functional form of the NLO formula (solid), which comes from 
the chiral logarithm originating from pion loops. The NLO formula works precisely up to 
A ~ 0.025, and the deviation is still within two sigma at A ~ 0.04, which is about the half 
of the physical strange quark mass mf^^^. This is a typical range where the NLO ChPT is 
valid, and the higher order corrections would become sizable above this value. In contrast 
to a recent work by Giusti and Liischer 36| , where they take a wider range of A (up to A ~ 
95 MeV) into the analysis, we conservatively choose the reference points below 30 MeV, 
so that (partially quenched) ChPT with an imaginary valence quark mass iX can be safely 
applied. 

In Figure [H the difference between the Nf = 3 and Nf = 2 (dashed curve) formulas is 
not sizable below A ~ 0.03-0.04. This is natural because the strange quark (with = 0.08) 
decouples from the dynamics of the low-lying modes. As a result, the extraction of Seg does 
not significantly depend on the formula we use [Nf = 3 or Nf = 2). 

The convergence of the chiral expansion is better in the e regime as shown in Figure [21 in 
which the Nf = 2 + 1 lattice data at rriud = 0.002 and rris = 0.080 are plotted. In the plot of 
the mode number (lower panel), the LO and NLO curves coincide up to A ~ 0.025. Beyond 
this value, we observe some deviation, which is also seen in the histogram plot (upper panel). 

The NLO ChPT correction in this work explains the disagreement of the lattice data 



with the expectation from the random matrix theory found in our previous work [15 



16|. 



Namely, if we adjust the parameter Seg using the lowest eigenvalue distribution (the first 
peak of the histogram), then the second peak would not agree at the leading order. Indeed, 
the NLO contribution is responsible for this. 
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A^/=3 ChPT (reduced) Nf=2 ChPT 
Nf (lattice) rriud ms comment 

Seff F Eeff F 



2 0.002 oo 




0.00218(19) 


0.059(65) 


(0= 


2.35) 


0.015 oo 




0.00362(15) 


0.0527(20) 






0.025 oo 




0.00353(15) 


0.0664(90) 










n (^(^'^R9(^ a) 

U.UUOO^l ±^ ) 


n Dfisi (fid') 






0.050 oo 




0.00449(15) 


0.0644(20) 






0.050 oo 




0.00400(16) 


0.0728(60) 


(0 


=-2) 


0.050 oo 




0.00482(19) 


0.0636(19) 


(0 


=-4) 


0.070 oo 




0.00480(15) 


0.0707(23) 






0.100 oo 




0.00478(12) 


0.0862(72) 






2+1 0.002 0.080 0.00204(07) 


0.0469(102) 0.00204(05) 


0.0425(49) 






0.015 0.080 0.00314(18) 


0.0536(15) 


0.00305(17) 


0.0551(16) 






0.015 0.080 0.00354(48) 


0.0521(25) 


0.00319(58) 


0.0558(62) 


(O 


= 1) 


0.015 0.080 0.00273(06) 


0.0520(25) 


0.00270(06) 


0.0545(26) 


iL-- 


=24) 


025 080 00333fl8) 


0.0624(20) 


0.00326(18) 


0647f20) 






0.025 0.080 0.00299(06) 


0.0600(23) 


0.00297(05) 


0.0629(24) 


(L-- 


=24) 


U.UOO U.UOU U.UU^U^lOc/l 


0.0636(17) 


0.00393(36) 


U.UUUUl -LU J 






0.050 0.080 0.00423(22) 


0.0696(16) 


0.00413(21) 


0.0738(16) 






0.015 0.100 0.00309(14) 


0.0564(19) 


0.00303(13) 


0.0578(19) 






0.025 0.100 0.00349(20) 


0.0622(17) 


0.00342(19) 


0.0642(17) 






035 100 00418f40') 


0.0647(14) 


0.00409(38) 


0673fl4~) 






0.050 0.100 0.00383(13) 


0.0713(16) 


0.00376(13) 


0.0747(16) 






3 0.025 0.025 0.00335(23) 


0.0531(10) 










0.035 0.035 0.00334(21) 


0.0612(24) 










0.080 0.080 0.00453(23) 


0.0767(14) 










0.100 0.100 0.00520(22) 


0.0835(22) 











TABLE II: Extracted values of and F using iV/ = 3 and iV/ = 2 CtiPT. 
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0.005 
0.004 




FIG. 1: Nf =2+1 lattice QCD results for the spectral density Pq(A) (top panel) and the mode 
number Nq{X) (bottom panel) of the Dirac operator at rriy^d = 0.015, = 0.080 and Q = 0. 
The lattice data (histogram (top) or solid symbols (bottom)) are compared with the NLO ChPT 
formula drawn by solid curves. For comparison, the prediction of the leading-order e expansion 
(dotted curves) and that of the NLO formula but with Nj = 2 flavors (dashed) are also shown. 
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FIG. 2: Same as Figured! but at rriud = 0.002 and nis = 0.080. 
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V. SCALING TESTS WITH VARIOUS PARAMETERS 



A. Sea quark masses 

Figures E] and m compare the spectral density and its integral at various sea quark masses. 
The data are obtained for Nf = 2 (Figure |3]) and Nf = 2 + 1 (Figure H]) with up and down 
quark masses in the p regime {rriud = 0.050, 0.035, 0.015) and in the e regime (m^^ = 0.002). 
Note that the /3 value of the e regime run in Nf = 2 QCD is slightly higher (/3 = 2.35) than 
in other runs (/3 = 2.30). In the plot, we adjust the value of A by a factor of 1.065 which 
corresponds to the ratio of lattice spacings between the two /3 values. 

In the plots, we clearly observe the sea quark mass dependence. For heavier sea quarks, 
the spectral density shows a higher peak near the lowest eigenvalue (around A ~ 0.004), and 
the peak height becomes lower by reducing the quark mass. As one enters the e regime, the 
lowest eigenvalue is pushed up to A ~ 0.015. This is what we expect for the effect of the 
fermionic determinant Jlfcl-^fc + ^tdf' ■ Namely, when the quark mass is reduced to the value 
around (or below) the lowest eigenvalue, those eigenvalues are suppressed. 

The expectation from NLO ChPT precisely follows the lattice data (solid curves in Fig- 
ures [3] and H]) . Here the values of Scg and F are determined for each set of sea quark mass, 
so that the comparison is not parameter-free. But, still the precise agreement of the shape 
of the spectral density is encouraging. We investigate the mass dependence of Sefr in the 
next section. 

In Figure similar plots are also shown for the Nf = 3 lattice data where three sea 
quarks have a degenerate mass (m„^ = m^). We have four values of the quark mass (0.100, 
0.080, 0.035 and 0.025) in the p regime. In Nf = 3 QCD, we observe a stronger dependence 
on the quark mass than in Nf = 2 or 2 -|- 1, as suggested by the NLO formula for Egg in 

m- 

B. Topological charge 

In the low eigenvalue region, the Dirac spectral density is known to be sensitive to the 
topological charge Q of the gauge fields, which is clearly seen in Figures and [3 The 
solid curves in the plots show the expectation from the NLO ChPT with input parameters 
determined from the Q = lattice data. Therefore, there is no free parameter to be adjusted 
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in the curves for the Q 7^ cases. We observe that the Q dependence of the lattice data is 
quahtatively well described by ChPT. 

Note, however, that the extracted values of Sefr and F from Q 7^ data show a 2.2 a 
difference for Nf = 2 while they are consistent in the _ Nj_ = 2 + 1 case. Since the topological 



charge dependence is a part of finite volume effects |37H39l|. which should be accounted for 
by the effective theory analysis, the deviation suggests possible higher order effects in the 
1/L expansion of ChPT. Indeed, the Nf = 2 runs are carried out with a shorter temporal 
extent T = 32 than that of Nf = 3 (T = 48), so that the lattice volume is ~ 17% smaller in 
the physical unit. Higher order finite volume effects may therefore be enhanced for Nf = 2. 
In the final results of the Nf = 2 data, we add this ~11% deviation as an estimate of the 
systematic error due to the finite volume. 



C. Finite volume 



The finite volume scaling can be tested more explicitly with the Nf = 2 + 1 runs, for which 
the L = 24 (~ 2.7 fm) lattice data are available. Here we note that the comparison of Eeff 
obtained on the L = 16 and L = 24 lattices is not straightforward, because there is a finite 
volume effect encoded in ^i(M^) in the definition of T^^s ( ITT|) . It is still possible to analytically 
convert the values of ^i(M^) in different volumes. The results for J^^s obtained on the L = 24 
lattice are converted to those at the smaller volume as Seg = 0.00273(06) — i- 0.00305(15) at 
mud = 0.015 and Ees = 0.00299(06) -> 0.00336(12) at mud = 0.025, which agrees with the 
values calculated on the L = 16 lattice: Sefr = 0.00314(18) and 0.00333(18), respectively. In 
fact, as shown in Figure [H we find that the same inputs of (converted) Ses and F describe 
the data at different volumes very well. 

From these analysis, the systematic error due to finite volume is estimated as ~ 3%, 
which is taken from the difference between the L = 16 and L = 24 (after the conversion) 
results. 
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0.01 0.02 0.03 0.04 0.05 0.06 




FIG. 3: Nf = 2 lattice QCD results for the spectral density Pq{X) (top panel) and the mode 
number Nq{X) (bottom panel) of the Dirac operator at various sea quark masses. The global 
topological charge is fixed to zero. The NLO ChPT (Nf = 2) results are drawn by solid curves. 
Note that the (3 = 2.35 data are rescaled as A — t- 1.065A and p — t- 1.209p according to the difference 
of the lattice spacing a. 
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FIG. 4: Same as Figure [3] but for the Nj = 2 + 1 data at = 0.080. Solid curves show the NLO 
ChPT (Nf = 3) formula. 
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m„^=0.056, Q= [ 
m„^=0.050, Q=-2 [ 
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FIG. 6: Dependence on the topological charge of the Nf = 2 lattice QCD results at rriud = 0.050. 
For the ChPT curves (solid), the same input values of Seg, F determined from Q = are used. 
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FIG. 7: Same as Figure E] but for A^/ = 2 + 1 QCD at ruud = 0.015 and = 0.080. 
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FIG. 8: Spectral density (top) and the mode number (bottom) on the larger volume lattice (L = 24 
[~2.7 fm]) at rriud = 0.015 and 0.025 with a fixed value of = 0.080. For comparison, L = 16 
lattice results as well as the ChPT predictions (solid curves) are shown in the bottom panel where 
the same input values of Sgfr, F (from L = 24 results) are used. 
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VI. DETERMINATION OF THE CHIRAL CONDENSATE 



The extracted values of T,cs and F for each lattice ensemble are summarized in Table ITTl 
Note that Seg is extracted at the NLO accuracy, while the value of F which first appears in 
the NLO term, might receive larger systematic corrections from the NNLO contributions. 

As already noted above, since rris is fixed at a large value (0.080 or 0.100) in the Nf = 2 + 1 
ensembles, there is little difference between the reduced Nf = 2 and Nf = 2 + 1 ChPT 
analysis near the chiral limit of rriud'- ^es and F are almost equal within the statistical 
errors. The difference between = 0.080 and rris = 0.100 is also small (always less than 
la), and therefore we concentrate on the data at rris = 0.080 in the following. 

Now we analyze the sea quark mass dependence of Seg. The NLO formula for T^^s ffTT]) 
contains the low energy constants S, F and Lq as parameters. The chiral condensate S, 
in particular, appears at the leading order, and its determination through the quark mass 
dependence of Seg is valid at the NLO accuracy, while the other parameters controlling the 
NLO correction can be determined at the leading-order. 

In the fit of the lattice data, we attempt two procedures: a 3-parameter (S, F, Lq) fit 
without any additional inputs, and 2-parameter (S, Lq) fits for several input values of F. 
In the 2-parameter fits, the value of F determined in our previous works both in the p 
regime (F = 0.0474(30)) 4^ and in the e regime (F = 0.0524(34)) 41| are used for the 
Nf = 2 lattice ensembles. Any difference due to the different input values of F suggests 
some systematic error (although it turns out to be only ~ l.lcr in the final result). For 
the Nf = 2 + 1 and 3 runs, we use a naive (linear) chiral limit of F given in Table [TTl of 
which values are F = 0.0411 for Nf = 3 ChPT and F = 0.0406 for reduced Nf = 2 ChPT, 
respectively. 

The chiral extrapolation of in Nf = 2 QCD is shown in Figure [91 From the plot, 
we can see the crucial role played by the data point at rriud = 0.002 which is in the e 
regime. Without this data point, one may naively expect that the data do not have enough 
sensitivity to probe the curvature due to the chiral logarithm and the chiral extrapolation 
favors a larger value of S (~ 0.0032). Taking the e regime data into account, the presence 
of chiral logarithm is consistent with the negative curvature seen in the data. 

The extracted values of the LECs in Nf = 2 ChPT are summarized in Table IIIII We 
attempt the 2-parameter fits of various number of data points, 3-7, taken from the lowest 
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FIG. 9: Chiral extrapolation of Sgfr in Nf = 2 QCD. The data point in the e regime {rriud = 
0.002) is rescaled to absorb the small difference of the lattice spacing. The two-parameter fits with 
an input F = 0.0474 |4(| for various number of data points included in the fit are drawn together 
with the lattice data points (open circles). The chiral limit is that for the 3-point fit. 

quark mass and the 3-parameter fits with 4-7 data points. In the table, the range of m^^ 
used in the fit is listed. For the 2-parameter fits, we take two input values of F. The quality 
of the fits can be inferred from the value of x^/d.o.f. also listed in the table. 

As far as the heaviest point is discarded in the fits, the resulting value of S is insensitive 
to the input value of F, and it is consistent with the three-parameter fit as well. On the other 
hand, the determination of Lg is unstable, but all the data suggest |Lg| < 0.001. We take 
S = 0.00246(15) and Lg = —0.00009(13) as the central values, which are from the 4-point 
fit with the input F = 0.0474. For the final results of this paper, we take the deviation from 
the other input value of F, as well as the other fitting ranges, as a systematic error due to 
the chiral extrapolation. 

For the Nf = 2 + 1 lattice data, the chiral extrapolation is more stable because of the 
precise data point in the e regime, which is simply due to higher statistics we accumulated. 
Figure [10] clearly shows the logarithmic curvature. Namely, a naive linear extrapolation of 
four data points in the p regime {rriud = 0.015-0.050) would lead to ~ 0.0028 in the chiral 
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mud fit range 


Nf= 

T, 


=2 ChPT LECs 

F LI 


xVd.o.f. 


2prm fit 


{F =0.0474) 








0.002-0.025 


0.00243(18) 


[0.0474] 


-0 00004f21) 


4.8 


0.002-0.035 


0.00246(15) 


[0.0474] 


-0.00009(13) 


2.5 


0.002-0.050 


0.00233(12) 


[n 04741 


OOOOTH 0^1 


2.4 


0.002-0.070 


0.00233(09) 


[0.0474] 


00006f07) 


1.8 


0.002-0.100 


0.00246(07) 


[0.0474] 


-0.00004(03) 


2.4 


2prm fit 


{F =0.0524) 








0.002-0.025 


0.00250(19) 


IQ 05241 


00012('30') 


5.2 


0.002-0.035 


0.00255(15) 


[0.0524] 


00002(18) 


2.7 


0.002-0.050 


0.00243(12) 


[0.0524] 


0.00021(14) 


2.4 


0.002-0.070 


0.00246(10) 


[0.0524] 


0.00018(09) 


1.8 


0.002-0.100 


0.00262(08) 


[0.0524] 


0.00001(04) 


2.8 


3prm fit 










0.002-0.035 


0.00174(45) 0.0290(71) 


-0.00024(03) 


2.5 


0.002-0.050 


0.00246(68) 


0.0547(78) 0.00038(90) 


3.5 


0.002-0.070 


0.00237(38) 0.0489(15) 0.00009(35) 


2.4 


0.002-0.100 


0.00206(26) 0.0386(48) 


-0.00009(02) 


2.3 



TABLE III: Nf = 2 lattice results for S, F and LQ{fisub = 770 MeV) extracted comparing with 
the Nf = 2 ChPT formula. The results for 2- and 3-parameter fits are listed. The values in the 
parenthesis [• • •] are used as an input of the chiral fit. 

limit, while the e regime point is lower (~ 0.0020). 

We analyze the Nf = 2 + 1 lattice data listed in Table |TT] using the 2- and 3-parameter 
fits. We consistently use the Nf = 3 formula for the data obtained with Nf = 3 ChPT 
(fourth column in Table HTl) . The same applies for the reduced Nf = 2 ChPT analysis. The 
fit results are presented in Table [IV] and in Table |Vl The curves in Figure [TO] represent the 
Nf = 3 fits for various numbers of data points. We find that all the curves go through the 
data points except for the heaviest one. The chiral limit shown by a cross symbol is almost 
unchanged by taking different fit schemes (2-parameter or 3-parameter) and the number of 
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FIG. 10: Same as Figure [9] but for Nj = 2+1 QCD. The fit curve corresponds to that of the 
3-parameter fit. 



Nf=3 ChPT LECs 
rriud fit range x^/d-o.f. 



2prm fit 








0.002-0.025 


0.00186(09) [0.0411] 


-0.00013(09) 


1.0 


0.002-0.035 


0.00186(09) [0.0411] 


-0.00014(08) 


0.6 


0.002-0.050 


0.00186(09) [0.0411] 


-0.00014(07) 


0.5 


0.002-0.080 


0.00185(08) [0.0411] 


-0.00014(07) 


0.4 


3prm fit 








0.002-0.035 


0.00185(10) 0.0433(13) 


-0.00023(53) 


1.3 


0.002-0.050 


0.00186(09) 0.0406(05) 


-0.00012(25) 


0.7 


0.002-0.080 


0.00186(08) 0.0413(02) 


-0.00015(09) 


0.5 



TABLE IV: Nf = 2 + 1 lattice results for S, F and Lq {^sub = 770 MeV) extracted using 
the A'^^ = 3 ChPT formula. The results for 2- and 3-parameter fits are listed. The values with 
parenthesis [• • •] are used as an input for the chiral fit. SP'^y^ denotes the rriud = and = oo 
hmit of Seff with mg = 0.080 fixed. 
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reduced Nf=2 ChPT LECs 

rriud fit range x^/d-o.f. 
— — - J. 

^6 



T, F LI 



9r)rTn fit, 








0.002-0.025 


0.00199(06) [0.0406] 


-0.00044(10) 


0.9 


0.002-0.035 


0.00197(06) [0.0406] 


-0.00040(09) 


0.7 


002-0 0,^0 


00197l'05'l Co 04061 


-0 0003qi'06'l 


4 


0.002-0.080 


0.00198(05) [0.0406] 


-0.00042(04) 


0.5 


3prm fit 








0.002-0.035 


0.00197(09) 0.0407(10) 


-0.00042(51) 


1.3 


0.002-0.050 


0.00198(07) 0.0416(05) 


-0.00038(21) 


0.6 


0.002-0.080 


0.00196(07) 0.0399(03) 


-0.00044(07) 


0.5 



TABLE V: Same as Table HV] but extracted using reduced Nf = 2 ChPT. 

data points included in the fit (4, 5, or 6). This is because the precise e regime point works 
as an anchor near the chiral Hmit. 

For Nf = 3 ChPT, 1]^^^' hsted in Table HV] denotes the value of in the limit of 
f^ud and y — 7- oo with a fixed strange quark mass = 0.080. This corresponds to the 
physical value of the chiral condensate —{uu) = —{dd) defined in the limit of vanishing up 
and down quark masses. 

From Tables IIVI and IVj we can see that results for S^'^^^ obtained via Nf = 2 + 1 and 
reduced Nf = 2 ChPT formulas are consistent with each other. (We assume that T,^^^^ in 
the A^j = 2 + 1 theory corresponds to E in reduced Nf = 2 ChPT.) The result is also stable 
against the changes of the number of fitting parameters and the fitting range. For the other 
parameters, F and Lg, we find larger dependence on the choice of fit procedures, which is 
expected because they only appear in the NLO terms. 

We take Sp^^^ = 0.00186(09), F = 0.0406(05) and = -0.00012(25) as the central 
values, which are obtained from the 5-point fit with three free parameters in Nf = 3 ChPT. 
As mentioned above, the other results are used to estimate systematic errors. 

For the degenerate {rriud = itLs) Nf = 3 lattice results, the number of data points is 
limited to four. The 3-parameter fit, therefore, does not work and we restrict ourselves to 
the 2-parameter fit. We attempt the fit with two values of F, 0.0431 and 0.0531, as the input. 
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FIG. 11: Same as Figure [9] but for the degenerate {rriud = rng) Nf = 3 QCD. Results of the 
2-parameter fits and a combined Nf = 3 fit are plotted. 



ruud, ms fit ran^ 


iV/=3 ChPT LECs 

s spi^y^ F 




xVd.o.f. 


2prm fit 








0.025-0.100 


0.00152(11) 0.00204(04) [0.0431] 


0.00010(10) 


3.1 


0.025-0.100 


0.00182(14) 0.00242(06) [0.0531] 


0.00031(17) 


2.8 


3prm fit 


(combined with Nf = 2+1 data) 






0.002-0.100 


0.00145(12) 0.00191(08) 0.0401(17) 


0.00003(7) 


1.7 



TABLE VI: The degenerate Nf = 3 lattice results for the LECs of iV/ = 3 ChPT. Here, S denotes 
the chiral condensate at the limit m^^ = = while T,^^^^ is the one with = 0.080 fixed. 

The former value is an Nf = 3 chiral limit of F in Table [IT] taken with a linear function in 
the quark mass and the latter is the one at the lightest sea quark mass rriud = f^s = 0.025. 
Due to the lack of the e regime data point, the chiral limit is not as stable as in the Nf = 2 
or Nf = 2 + 1 data. In fact, the resulting value of S strongly depends on the input value 
of F. Between the two representative values of F, S changes about 20%. Since F controls 




combined fit (w/ Njr= 2+1 lattice) 

(chiral limit) i- -x- - - 

_i I I I I 



30 



the NLO effects as seen in (fTTI) . this change suggests that the one-loop calculation is not 
sufficient to control the rriud = dependence. 

We also attempt a combined fit of all the Nf = 2 + 1 and Nf = 3 data points using 
the Nf = 3 ChPT formula. The total nine data points are simultaneously fitted with a 
reasonable x^/d.o.f. 1.7). The result is given in Table IVT] and plotted in Figure [TT] (and 
[TOj) which indicate a strong effect of chiral logarithm in the Nf = 3 data compared to the 
case of Nf = 2 shown in Figure |9l The chiral limit is less than half of the value at the lowest 
p regime point {rriud = '"^s = 0.025). One should note, however, that the data points in the 
fit includes those with = 0.080 and 0.100 which are out of the typical convergence region 
of chiral expansion < so that the result may contain large systematic effects. 

It is still remarkable that all the results suggest that the chiral condensate S of Nf = 3 
QCD is smaller than Sp^^" of A^/ = 2 + 1 or S of A^/ = 2 QCD. This is consistent with the 
view that the chiral condensate decreases and eventually disappears as the number of flavor 
increases and the asymptotic freedom is lost. We take S = 0.00145(12), F = 0.0401(17) and 
Lg = 0.00003(7) determined from the combined fit as the central values of Nf = 3 ChPT 
parameter, and will include a ~ 26% deviation from the 2-parameter fit with F = 0.0531 as 
a systematic error in the final results. 

So far we have treated the strange quark mass fixed at = 0.080 in the Nf = 2 + 1 
studies. In fact, the combined fit of the degenerate Nf = 3 and 2+1 results implies that 
Soff in the chiral limit of rriud changes by less than 1% when rris varies between 0.060 and 
0.100. We can, therefore, safely ignore the error due to a slight mismatch of the strange 
quark mass from its physical value. This weak sensitivity to supports the assumption 
that the strange quark at the physical mass is almost decoupled from the low energy theory 
and the use of the reduced Nf = 2 ChPT formula to fit the Nf = 2 + 1 lattice QCD data. 

VII. SUMMARY AND CONCLUSION 

Before quoting the final results, let us discuss the possible systematic errors. 

Our simulation uses an exactly chiral-symmetric Dirac operator and has reached almost 
the chiral limit: 771^4^=0.002 for the Nf = 2 and 2+1 runs, which corresponds to ~ 3 MeV 
in the physical unit. The NLO ChPT formula is valid in both the e and p regimes. As 
a result, the chiral extrapolation of Sgfr in Nf = 2 and 2+1 QCD is stable. As discussed 
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in the previous section, by varying the fit range (and Nf in the ChPT formula for the 
Nf = 2 + 1 analysis), we estimate the systematic effects due to the chiral extrapolation as 
(+3.7%/- 5.2%) for S of A^/ = 2 QCD and (+6.5%/- 0.5%) for Sp^^^^ of A^/ = 2 + 1 QCD, 
respectively. The upper and lower limits come from the variation of S with the various fit 
schemes. We also found that the dependence is negligible in the range = 0.060-0.100, 
so that SP*^^*^ can be treated as the chiral condensate at the physical value of the strange 
quark mass. 

On the other hand, for the degenerate Nf = 3 lattice data, our estimate for the systematic 
error in the chiral fit is larger because of the bad convergence of Nj = 3 ChPT and smaller 
number of data points. In fact, we observe that S moves as large as 26% depending on the 
fit schemes, from which we estimate the systematic error from this source to be ±26% for 
S in the Nf = 3 chiral limit. 

The finite volume effects have also been discussed in the previous sections. For the Nf = 2 
lattice results, we expect that a possible higher order effect in the 1/L expansion beyond 
one-loop ChPT is partly reflected in the difference among different topological sectors. We 
thus estimate the systematic error from this source to be ±11%. For the Nf = 2 + 1 (and 3) 
case, we use more direct comparison with L = 24 lattice results and the systematic error is 
estimated as (+0.9%/— 2.9%). With a naive order counting, the leading finite volume effect 
is estimated as 1/{F'^V), which is the two- loop effect in the p expansion. With F = 71 MeV 
(see below), this gives a large value (~ 0.52) as the size of the two-loop correction. In fact, 
including the numerical coefficient /3i given in the Appendix |A] the one-loop correction is 
/3i/(-F^V"^/^) ~ —0.06. If we assume that the numerical coefficient is also small (~ 0.05) at 
the two-loop order, this naive estimate gives a 3% effect, which is in the same ball park as 
the estimate given above. The small numerical coefficient at the two-loop level is indeed 
obtained in a recent study [42]. 

Since our lattice studies are done at only one value of (3, it is difficult to estimate the 
size of discretization effects. It should be partly reflected in the mismatch of the observables 
measured in different ways. For instance, the inverse lattice spacing determined from the 
0-meson mass, 1/a = 1.774(17) GeV, is 1% larger than the determination from the fi-baryon 
mass, 1/a = 1.759(10) GeV 35|. The latter corresponds to the Sommer scale tq = 0.51 fm. 



which is higher than the nominal value 0.49 fm or the recently favaored value 0.46-0.47 fm by 
about 4-10%. On the other hand, a naive order counting (oAqcd)^ with Aqcd ~ 450 MeV 
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Nf = 2 ChPT = 3 ChPT 







pphys 


i/ d 








renormalization 


+ 1.4 0/ 
-1.1 /o 


+1.2 

-1.1 


% 


+ 1.2 07 
-1.1 /o 






chiral fit 


+ 1.2 0/ 


+2.2 
-0.2 


% 


±8.7 % 


±8.0 


% 


finite volume 


±3.7% 


+0.3 
-1.0 


% 


+0.3% 
-1.0 


±3.0 


% 


finite a 


±7.0 % 


±7.0 


% 


±7.0 % 


±7.0 


% 


total 


+8.1 0/ 
-8.2 /o 


+7.4 
-7.2 


% 


±11 % 


±11 


% 



TABLE VII: Systematic errors for [S^S(2 GeV)]i/3 and F. The total errors are obtained by adding 
each estimate by quadrature. 

suggests a systematic effect of ~ 7%, which is consistent with the above mismatch. We 
therefore add this naive estimate, ±7%, as the systematic error due to finite lattice spacing. 

We convert the value of the condensate to the definition in the standard renormaliza- 
tion scheme, i.e. the MS scheme. By using the nonperturbative renormalization tech- 



nique through the RI/MOM scheme we obtained the Z factor in our previous work 43|] as 
l/Zsi2GeV) = 0.804(10) (1^1) for Nj = 2 and 0.792(10) (1^^) ioi Nj = 2 + 1 and 3, where 
the errors are statistical and systematic, respectively^. 

Including all the systematic effects above, which are summarized in Table IVIIt we obtain 
the chiral condensate of up and down quarks in their massless limit as 



E^'^(2 GeV) 



[242(05) (20) MeV]=^ for Nf = 2 QCD 

[234(04) (17) MeV]3 for Nj = 2+1 QCD (at physical m,) , (15) 
^ [214(06) (24) MeV]3 for A^^ = 3 QCD (at m, = 0) 

where the errors are again statistical and systematic, respectively. Here, the total systematic 
errors are obtained by adding each estimate by quadrature. Note that the result for Nf = 
2 ± 1 is slightly changed from |18j because a different input for the scale determination is 
used. We find a nontrivial dependence on the chiral condensate: as the strange quark 
mass goes down from rris = oo {Nf = 2 QCD) to the chiral limit rris = 0, the value of E 
decreases. 



2 



The value for Nf = 2 + 1 is slightly changed from [43[ due to the different determination of the lattice 



scale, that affects the renormalization group running of the Z factor. 
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The chiral condensate determined in this work ( !T5|) is consistent with those in our pre 



vious results obtained from the pseudoscalar meson mass 



20 



regime ChPT 



40| and from the topological 



susceptibility j44446|. A similar work done using the Nf = 2 Wilson fermion and the p 



36j quoted J:^^{2 GeV) = 276(3) (4) (5) MeV, which is slightly higher than 



our result. A recent work [47| in a mixed action approach (overlap valence + Wilson sea) 
has also reported a larger value. More detailed study would be necessary to understand the 
source of the discrepancy, if it is significant. 
From the NLO terms, we also obtain 

F = 71 (3) (8) MeV for A^^ = 3 QCD, (16) 

(or ^/2F =100(4) (11) MeV) and 

f -0.00009(13) (30) for A^. = 2 QCD, ^ ^ 

L^(770MeV) = <^ v A ; / ^ ^ ^^^^ 

[ 0.00003(07) (17) for Nf = 3 QCD, 

where the systematic errors are estimated in a similar manner. For F, their estimates are 
listed in Table lylTl while for Lg, the systematic error is dominated by the one from chiral 
extrapolation, as seen in Table UTTl and II Vl Although the accuracy for these quantities is not 
as good as that of the chiral condensate, they provide important consistency checks. 

In this study, we have investigated the eigenvalue spectrum of the QCD Dirac operator, 
which is free from ultraviolet power divergences. The lattice QCD results show a good 
agreement with the ChPT calculation at NLO in the region of A less than mfy''/2. In 
particular, the effect of pion-loop, or the chiral logarithm, is clearly seen. The dependence 
on the volume V, the quark masses rriud and m^, and the topological charge Q is also well 
described by ChPT. Result for the chiral condensate extracted from this study is therefore 
robust, as the systematic errors are controlled except for that coming from the discretization 
effect. 

Our work has also addressed a nontrivial flavor dependence of the chiral condensate. As 
the strange quark mass is reduced towards the chiral limit, its dynamical effect is seen as 
lowering the value of S. 
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Appendix A: Shape coefficient /3„ 

In ( ITOj) . we need to calculate the coefficients /3„'s jsol, which depend only on the shape 
of the four-dimensional box. The definition of /3„ is given by 

2 \ , , , ^ as -ln47r + 7-3/2 , 



Yl/2\3 /yl/2^ 



oo 

S{x) = J2 expi-irk^x), (A3) 

fc=— oo 
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T/L /3i /32 ^3 ^4 ^5 /36 

1 0.1405 -2.030x10-2 -4.820x10^^ 2.531x10-^ -2.238x10"*^ 2.672x10^'^ 

2 0.08360 -1.295x10-2 -1.778x10-^ 3.265x10"'' -9.120x10"^ 3.250x10"^ 

3 -0.04194 0.01215 -9.508x10-^ 3.622x10-^ -1. 898x10^3 1.248x10^3 

TABLE VIII: Numerical results for /3„ for T/L = 1, 2 and 3. 

where 7 ~ 0.577215665 ■■ ■ is the Euler's constant. Here the summation in S{x) is well 
approximated by a truncation \k\ < 20. For the case with T/L =1,2 and 3, the numerical 
values of /3„ are listed in Table IVIIII 
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